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TOWARDS THE MODELING OF NEURONAL FIRING BY GAUSSIAN 

PROCESSES 

E. Di Nardo, A.G. Nobile, E. Pirozzi and L.M. Ricciardi 



Abstract. This paper focuses on the outline of some computational methods for 
the approximate solution of the integral equations for the neuronal firing probabil- 
ity density and an algorithm for the generation of sample-paths in order to construct 
histograms estimating the firing densities. Our results originate from the study of 
non-Markov stationary Gaussian neuronal models with the aim to determine the neu- 
ron's firing probability density function. A parallel algorithm has been implemented 
in order to simulate large numbers of sample paths of Gaussian processes character- 
ized by damped oscillatory covariances in the presence of time dependent boundaries. 
The analysis based on the simulation procedure provides an alternative research tool 
when closed-form results or analytic evaluation of the neuronal firing densities are not 
available. 



1 Introduction 

This contribution deals with the implementation of procedures and methods, worked out 
in our group during the last few years, in order to provide algorithmic solutions to the 
problem of determining the first passage time (FPT) probability densities (pdf) and its 
relevant statistics for continuous state-space and continuous parameter Gaussian processes 
describing the stochastic modeling of a single neuron's activity. 

In most modeling approaches, it is customary to assume that a neuron is subject to 
input pulses occurring randomly in time, (see, for instance, \\?>\ and references therein). As 
a consequence of the received stimulations, it reacts by producing a response that consists 
of a spike train. The reproduction of the statistical features of such spike trains has been 
the goal of many researches who have focused the attention on the analysis of the interspike 
intervals. Indeed, the importance of the interspike intervals is due to the generally accepted 
hypothesis that the information transferred within the nervous system is usually encoded 
by the timing of occurrence of neuronal spikes. 

To describe the dynamics of the neuronal firing we consider a stochastic process X{t) 
representing the change in the neuron membrane potential between two consecutive spikes 
(cf., for instance, In this context, the threshold voltage is viewed as a deterministic 
function S(t) and the instant when the membrane potential crosses S(t) as a FPT random 
variable. 

The modeling of a single neuron's activity by means of a stochastic process has been the 
object of numerous investigations during the last four decades. A milestone contribution 
in this direction is the much celebrated paper by Gerstein and Mandelbrot [7] in which a 
random walk and its continuous diffusion limit (the Wiener process) was proposed with the 
aim of describing a possible, highly schematized, spike generation mechanism. However, 
despite the excellent fitting of a variety of data, this model has been the target of severe 
criticism on the base of its extreme idealization in contrast with some electrophysiological 
evidence: for example, this model does not take into account the spontaneous exponential 
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decay of the neuron membrane potential. An improved model is the so called Ornstein- 
Uhlenbeck (OU) model, that embodies the presence of such exponential decay. However, 
the OU model does not allow to obtain any closed form expression for the firing pdf, except 
for some very particular cases of no interest within the neuronal modeling context. Rather 
cumbersome computations are thus required to obtain evaluations of the statistics of the 
firing time. Successively, alternative neuronal models have been proposed, that include 
more physiologically features. The literature on this subject is too vast to be recalled here. 
We limit ourselves to mentioning that a review of most significant neuronal models can be 
found in ^H], Q3| and in the references therein. In particular, in jJJJ] it is presented an 
outline of appropriate mathematical techniques by which to approach the FPT problem in 
the neuronal context. 

We shall now formally define the firing pdf for a model based on a stochastic process 
X(t) with continuous sample paths. First, assume P[X(to) — xo] — 1, with xq < S(to), i.e. 
we view the sample paths of X(t) as originating at a preassigned state xq at initial time to. 
Then, 

T X0 = mf{t:X(t) >S(t)\, x < S(t Q ) 
t>t 

is the FPT of X(t) through S(t), and 

g[S(t),t\x M = §- t P(T X0 <t) 

is its pdf. 

Henceforth, the FPT pdf g[S(t), t\xo, to] will be identified with the firing pdf of a neuron 
whose membrane potential is modeled by X{t) and whose firing threshold is S(t). 

Throughout this paper, we shall focus our attention on neuronal models rooted on 
diffusion and Gaussian processes, partially motivated by the generally accepted hypothesis 
that in numerous instances the neuronal firing is caused by the superposition of a very 
large number of synaptic input pulses which is suggestive of the generation of Gaussian 
distributions by virtue of some sort of central limit theorems. 

It must be explicitly pointed out that models based on diffusion processes are character- 
ized by the "lack of memory" as a consequence of the underlying Markov property. However, 
in the realistic presence of correlated input stimulations, the Markov assumption breaks 
down and one faces the problem of considering more general stochastic models, for which 
the literature on FPT problem appears scarce and fragmentary. Simulation procedures 
then provide possible alternative investigation tools especially if they can be implemented 
on parallel computers, (see The goal of a typical simulation procedure is to sample N 
values of the FPT by a suitable construction of N time-discrete sample paths of the process 
and then to record the instants when such sample paths first cross the boundary. In such 
a way, one is led to obtain estimates of the firing pdf and of its statistics, that may be 
implemented for data fitting purposes. 

The aim of this paper is to outline numerical and theoretical methods to characterize the 
FPT pdf for Gaussian processes. Attention will be focused on Markov models in Section |2 
and on non-Markov models in Section[3] Finally, Section0]will be devoted to the description 
of some computational results. 

2 Gauss-Markov processes 

We start briefly reviewing some essential properties of Gauss-Markov processes. Let {X(t), t £ 
/}, where / is a continuous parameter set, be a real continuous Gauss-Markov process with 
the following properties (cf. [Hj): 

(i) mif) :— E[X(t)} is continuous in /; 
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(ii) the covariance c(s, t) :— E{[X(s) — m(s)] [X(t) — m(t)]} is continuous in I x I: 

(Hi) X(t) is non-singular, except possibly at the end points of / where it could be equal to 
to (i) with probability one. 



A Gaussian process is Markov if and only if its covariance satisfies 
(1) 



c(s,t)c(t,u) 
c(s, u) = TTT\ ' s ' t,U £ l,s < t < u. 



c(t,t) 

It is well known jS], that well-behaved solutions of are of the form 
(2) c(s,t) = h 1 (s)h 2 (t), s<t, 

where 

hi(t) 



(3) 



it) := 



h 2 (t) 



is a monotonically increasing function by virtue of the Cauchy-Schwarz inequality and 
hi(t) fi2(t) > because of the assumed nonsingularity of the process on /. The conditional 
pdf f(x, t\xo, to) of X(t) is a normal density characterized respectively by conditional mean 
and variance 



M(t\t ) = m(t) + 



h 2 {t) 
h 2 (t ) 



[x -m(t )} 



V(t\t 



h 2 (t) 



h x {t)-^-h x {t Q ) 



with t, to e/, to < t. It satisfies the Fokker-Planck equation and the associated initial 
condition 



df(x,t\y,r) 
dt 



[A 1 (x,t)f{x,t\y,T)] + \%- [A 2 {t) f(x,t\y,T)], 



d.r 



2 dx 2 



lim f(x,t\y,r) = 5(x - y), 

Tit 

with Ai(x,t) and A 2 (t) given by 

Ai(x,t) = m'(t) + [x- m(t)] 



W) 
h 2 (t) ' 



A 2 {t) = hl(t) r'(t), 



the prime denoting derivative with respect to the argument. 

The class of the Gauss-Markov processes {X(t),t 6 [0, +oo)}, such that f(x, t\y, r) 
f(x, t — r\y), is characterized by means and covariances of the following two forms: 

m(t) — Pit + c, c(s, t) = a 2 s + ci 

(0<s<t< +oo, /3i, c e R, ci > 0, a ^ 0) 



m{t) 



+ ce ftf , c(s,t) = ae 



thi 



c 2 e 



02 S 



2ci/3 2 



-/3 2 s 



0<s<t <+oo, /3i,c,c 2 £ R, <t ^ 0, ci ^ 0,/3 2 # 0, cic 2 - -5- > 
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The first type includes the Wiener process, used in [7j to describe the neuronal firing, while 
the second type includes the Ornstein-Uhlenbeck process that has often been invoked as a 
model for neuronal activity (see, for instance, 

Any Gaussian process with covariance as in (j2J can be represented in terms of the 
standard Wiener process {W(t),t > 0} as 

(4) X(t)=m(t)+h 2 (t)W[r(t)], 

and is therefore Markov. This last equation suggests the way to construct the FPT pdf of a 
Gauss-Markov process X(t) in terms of the FPT pdf of the standard Wiener process W(t). 
As an example, from (0J for the conditioned FPT pdf one has: 

(5) g[S(t),t\xo,t ] = ^ g w {S*[r(t)],r(t)\x* ,r(t )}, 

where r(t) is defined in J3Jl and gw[S*{i)), &\xq, i9 ] is the FPT pdf of W(-d) from Xg at time 
$0 to the continuous boundary with 

x -m[r-\& Q )\ q * ( „ v _ S[r- 1 {-&)]-m[r- 1 {d)) 

Results on the FPT pdf for the standard Wiener process can thus in principle be used 
via © to obtain the FPT pdf of any continuous Gauss-Markov process. For instance, if 
S*({>) is linear in gwlS*^), t9|xg, f?o] is known and g[S(t), t\xo, to] can be obtained via 
Instead, if gw ^o] is not known, a numerical algorithm, or a simulation 
procedure, should be used for the standard Wiener process and, after that, g[S(t), t\xo, to] 
can be obtained via the indicated transformation. 

The above procedure often exhibits the serious drawback of ensuing unacceptable time 
dilations (see 0]). For instance, exponentially large times are involved when transforming 
the Ornstein-Uhlenbeck process to the Wiener process, which makes such a method hardly 
viable. Hence, it is desirable to dispose of a direct and efficient computational method to 
obtain evaluation of the firing pdf. Along such a direction, in 0] it has been proved that 
the conditioned FPT density of a Gauss-Markov process can be obtained by solving the 
non-singular Volterra second kind integral equation 

g[S(t),t\x ,h] = -W[S(t),t\x ,to]+2 I g[S(r),T\x ,t ]^[S(t),t\S(T),T] dr 



to 



(6) (x < S(t )) 

with S{t),m(t),h 1 (t),h 2 {t) E C^I) and 

S'{t)-m'{t) S(t)-m(t) h[(t)h 2 (T) - h' 2 (t)h!(T) 



*[S(t),t\y, 



2 2 hx{t)h 2 {T) - h 2 (t)hi(r) 



l7 s y-m(r) h'^h^t) - h 2 (t)h[(t) \ 

(?) " h l{ t)h 2 {T)~h 2{ t)h l{ r) ) W®'* 1 ^ T] 

where f[x,t\y,r] is the transition pdf of X(t). Closed form solutions of arc available in 
4 for different families of boundaries. 

By making use of this result, in 0] an efficient numerical procedure based on a repeated 
Simpson's rule has been proposed to evaluate FPT densities of Gauss-Markov processes, 
that can be implemented to obtain reliable evaluations of firing densities for neuronal models 
based on Gauss-Markov processes. 
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3 Gauss non-Markov processes 

The methods proposed in the previous Section rest on the strong Markov assumption on the 
stochastic process modeling the neuron's membrane potential, which grants the possibility of 
making use of the mentioned analytic and computational methods for FPT pdf evaluations. 
This is not the case when the stochastic process used to model the neuron's firing mechanism 
is non-Markov. Here we shall focus our attention on Gauss non-Markov processes. However, 
thus doing we face the lack of effective analytical methods for obtaining manageable closed- 
form expressions for the FPT pdf, although some preliminary analytical results have been 
obtained by Ricciardi and Sato in for a class of stationary Gaussian processes. 

Indeed, if X(t) is one-dimensional non-singular stationary Gaussian process mean square 
differentiable, a series expansion for the FPT pdf was derived (see, In the following, 

for convenience, we shall take to = as initial time and, without loss of generality, assume 
that E[X(t)} = and P[X(to) — xq] — 1, with xq an arbitrarily specified initial state. 
Furthermore, the covariance function E[X(t)X(r)] :— j(t — r) will be assumed to be such 
that 7(0) = 1,7(0) = and 7(0) < (this last assumptions being equivalent to the mean 
square differentiable property). The FPT pdf of X{t) through S{t) is then given by the 
following expression 



(8) g[S(t),t\x } = W 1 (t\x )+J2(-l) i / dtj dt 2 --- dt z W l+1 

l=1 JO Jt x Jti„i 



(ti, . . .,ti,t\x ), 



with 

W„(ti, . . . ,t n \x ) 

pOO pOO n 

= I dzi--- dz n T[ [zi - S(ti)]p 2n [S(h), . . . , S(t n ); z 1 , . . . , z„\x ], 
JS(ti) Js(t n ) l=1 

where p 2 n{xi, . . . , x n ; Z\,. . . , z n \x ) is the joint pdf of {X(ti), . . .,X(t n ), Z{t x ) = X{t\), . . . , 
Z{t n )— X(t n )} conditional upon X(0) = xq. Due to the great complexity of the involved 
multiple integrals, expression (JHJ does not appear to be manageable for practical uses, even 
though it has recently been shown that it allows to obtain some interesting asymptotic 
results Since (JSJ) is a Leibnitz series for each fixed t > 0, estimates of the FPT pdf can 
in principle be obtained since its partial sum of order n provides a lower or an upper bound 
to g depending on whether n is even or odd. However, also the evaluation of such partial 
sums is extremely cumbersome. 

In conclusion, at the present time for this class of Gaussian processes, no effective 
analytical methods, nor viable numerical algorithms are available to evaluate the FPT 
pdf. A simulation procedure seems to be the only residual way of approach. 

To this aim, we have restored and updated an algorithm due to Franklin |F| in order to 
construct sample paths of a stationary Gaussian process with spectral density of a rational 
type and deterministic starting point. The idea is the following. Let us consider the linear 
filter 

p OO 

(9) X(t) = / h(s)W(t-s) ds 

Jo 

where h(t) is the impulse response function and W(t) is the input signal. By Fourier 
transformation, © yields 



(10) 



T x (u) = \H(u)\ 2 T w {lj) 
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Figure 1: Plot of the boundary S(t) given in for P = 0.5 and d = 0.25, 0.50 (bottom to top) 



where Tw{u) and Tx(<jj) are the spectral densities of input W(t) and output X(t), respec- 
tively, and where H(uj) denotes the Fourier transform of h(t). Equation itl"U|) is suggestive 
of a method to construct a Gaussian process X(t) having a preassigned spectral density 
rx(w) = r(a;). It is indeed sufficient to consider a white noise A(t), having spectral density 
Fiyfu) = 1, as the input signal and then select h(t) in such a way that \H(lu)\ 2 = T(ui). If 
the spectral density of X(t) is of rational type, namely if 



(11) 



7(*) e" 



dt 



P(iu) 



Q(iuj) 



where P and Q are polynomials with deg(P) < deg{Q), setting H(u>) — P(iu))/Q(iu), from 
fTTijl it follows 

where D = d/dt. To calculate the output signal X(t) it is thus necessary to solve first 
the differential equation Q(D) cj)(t) = A(t) to obtain the stationary solution <fi(t), and then 
evaluate X(t) — P(D) 4>{t). The simulation procedure is designed to construct sample paths 
of the process X(t) at the instants t = 0, At, 2 At, . . . where At is a positive constant time 
increment. The underlying idea can be applied to any Gaussian process having spectral 
densities of a rational type and, since the sample paths of the simulated process are gen- 
erated independently of one another, this simulation procedure is particularly suited for 
implementation on supercomputers. 

Extensive computations have been performed on parallel computers to explore the dif- 
ferent shapes of the FPT pdf as induced by the oscillatory behaviors of covariances and 
thresholds (cf., for instance, and |2])- 



4 Numerical comparisons 

The aim of this Section is to compare the behavior of the FPT pdf 's among Gauss-Markov 
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Figure 2: Plots refer to FPT pdf g(t) through the boundary IT51 with /3 = 0.5 and d = 0.25 for a 
zero-mean Gaussian process characterized by correlation function < 12ft . In Figure 2(a) given 
by l|14|l has been plotted. The estimated FPT pdf g(t) with a — 10 -10 is shown in Figure 2(b), 
with a = 0.25 in Figure 2(c) and with a = 0.5 in Figure 2(d). 

processes and Gaussian non-Markov processes, in order to analyze how the lack of memory 
affects the shape of the density, also with reference to the specified type of correlation 
function. For simplicity, set xo — and P[X(0) = 0] = 1, so that in the following we shall 
consider the FPT pdf g(t) := g[S(t),t\0, 0]. 

To be specific, we consider a stationary Gaussian process X(t) with zero mean and 
correlation function 

(12) 7 (t) := e"' 31 ' 1 cos(at), [3 e R + 

which is the simplest type of correlation having a concrete engineering significance 14 . 
When the correlation function is of type Q12[l . X(t) is not mean square differentiable, since 
7(0) 7^ 0. Thus the series expansion © does not hold. However, specific assumptions on 
the parameter a help us characterize the shape of the FPT pdf. 

We start assuming a = 0, so that the correlation function (|12fl can be factorized as 

7 (t) = e -P-r e -Ht-r) /3eR+,o<r<<. 

Hence, choosing h\(t) = e' 3 * and h,2{t) = e _/3t in J2J), X(t) becomes Gauss-Markov. There- 
fore, for any boundary S(t), the FPT pdf g(t) can be numerically evaluated by solving the 
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Figure 3: Same as in Figura^with d = 0.5. 



integral equation © . In the following, we consider boundaries of the form 



(13) S(t)=de-f >t h--^-ln 



i + iJl + 8exp 



4d 2 



s 2/3t _ I 



with rf > 0. It is evident that lim t ^o S(t) — d and that S(t) tends to as t increases. 
Furthermore, as d decreases, the boundary becomes flatter. In Figure ^ the boundary S(t) 
given in 113fl is plotted for (3 = 0.5 and for two choices of the parameter d, i.e. d = 0.25 
and d = 0.5. 

As proved in 0] for boundaries of the form l|13|) . the FPT pdf g(t) of a Gauss-Markov 
process admits the following closed form: 



(14) 



g(t) 



4d(3e 

- 1 



exp 



Ad 2 



3 2f3t _ l 



1 + 8 exp 



Ad 2 



3 2f3t _ l 



/[S(t),t|0,0] 



where f[S(t), t\0, 0] is the transition pdf of the Gauss-Markov process X(t). 

For a zero-mean Gauss-Markov process characterized by the correlation function Ijl2(l 
with (3 = 0.5 and a = 0, the FPT pdf g(t) given by through the boundary {T^J is 
plotted in Figure 2(a) for d = 0.25 and in Figure 3(a) for d = 0.5. Note that as d increases 
the mode increase, whereas the corresponding ordinate decrease. 
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Setting a ^ in 112[) . the Gaussian process X(t) is no longer Markov and its spectral 
density is given by 

d5) i»- m<s+<* + 0>) 



uj 4 + 2lu 2 {/3 2 - a 2 ) + {p 2 + a 2 ) 2 ' 

thus being of a rational type. Since in iJTSJl the de gree of the numerator is less than the 
degree of the denominator, it is possible to apply the simulation algorithm described in 
Section 3 in order to estimate the FPT pdf g(t) of the process. 

The simulation procedure has been implemented by a parallel FORTRAN 90 code on a 
128-processor IBM SP4 supercomputer, based on MPI language for parallel processing. The 
number of simulated sample paths has been set equal to 10 7 . The estimated FPT pdf g(t) 
through the boundary (|T3|> with /3 = 0.5 and d = 0.25 are plotted in Figures 2(b)-^2(d) for 
Gaussian processes with correlation function l|12(l having a = 10 -10 , 0.25, 0.5, respectively. 
For the same processes, Figures 3(b)-^3(d) show the estimated FPT pdf g(t) through the 
boundary l|13|l with j3 — 0.5 and d = 0.5. Note that as a increases, the shape of the FPT 
pdf g(t) becomes natter and the related mode increases. Furthermore, as Figures 2(a)-2(b) 
and Figures 3(a)-3(b) show, g(t) is very similar to g(t) for small values of a. 
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